1. Field of the Invention
This invention relates generally to the field of geophysical prospecting. More particularly, the invention relates to the field of seismic data processing. Specifically, the invention is a method for attenuating multiples in seismic data using a high-resolution Radon transform.
2. Description of the Related Art
In the field of geophysical prospecting, the knowledge of the subsurface structure of the earth is useful for finding and extracting valuable mineral resources, such as oil and natural gas. A well-known tool of geophysical prospecting is a seismic survey. A seismic survey transmits acoustic waves emitted from appropriate energy sources into the earth and collects the reflected signals using an array of receivers. Then seismic data processing techniques are applied to the collected data to estimate the subsurface structure.
In a seismic survey, the seismic signal is generated by injecting an acoustic signal from on or near the earth""s surface, which then travels downwardly into the earth""s subsurface. The acoustic signal may also travel downwardly through a body of water, in a marine survey. Appropriate energy sources may include explosives or vibrators on land and air guns or marine vibrators in water regimes. When the acoustic signal encounters a seismic reflector, an interface between two subsurface strata having different acoustic impedances, a portion of the acoustic signal is reflected back to the surface, where the reflected energy is detected by a receiver. Appropriate detectors may include particle motion detectors (such as geophones) on land and pressure detectors (such as hydrophones) in water regimes. Both sources and receivers may be deployed by themselves or, more commonly, in arrays.
The seismic energy recorded by each source and receiver pair during the data acquisition stage is known as a seismic trace. Seismic data traces contain the desired seismic reflections, known as the primary reflections or primaries. A primary reflection comes from the detection of an acoustic signal that travels from a source to a receiver with but a single reflection from a subsurface seismic reflector. Unfortunately, the seismic traces often contain many unwanted additional reflections known as multiple reflections or multiples, which can obscure and even overwhelm the sought-after primary reflections. A multiple reflection comes from the recording of an acoustic signal that has reflected more than once before being detected by a receiver. The additional multiple reflections could come from subsurface reflectors or from the surface of the earth in a land seismic survey and the water-earth or air-water interfaces in a water seismic survey. The recorded signals from multiples obscure the recorded signals from the primaries, making it harder to identify and interpret the desired primaries. Thus, the removal, or at least attenuation, of multiples is a desired step in seismic data processing in many environments. This is particularly so in marine seismic surveys, where multiples are especially strong relative to primaries. This is because the water-earth and, particularly, the air-water interfaces are strong seismic reflectors due to their high acoustic impedance contrasts.
One of many techniques well known in the art of seismic data processing for attenuating multiples is the use of a Radon transform. First, a single Fourier transform takes data from the space-time domain, or (x, t) domain, to the Fourier domain, or (x, xcfx89) domain. Here, co is the frequency. Then, for each frequency c, the Radon transform, or decomposition, transforms data to the Radon domain, or (q, xcfx84) domain. Here, q is a parameter associated with the curvature of the curvilinear trajectory describing the seismic event in the (x, t) domain and X (tau) is the zero offset intercept time. The Radon transform is given in simplest terms by the equation
Lm=d,xe2x80x83xe2x80x83(1) 
where L is a matrix representing the Radon transform, m is a vector representing the model being solved for in the Radon domain, and d is a vector containing the seismic data set being processed. The vector d is size M, where M is the number of available seismic traces di in the seismic data set d. Typically, M is the number of digitized seismic traces in a gather. The vector m is size N, where N is the number of spectral curvature components mi that d is being decomposed into by the Radon transform L. Thus, the matrix L representing the Radon transform is size Mxc3x97N and decomposes the seismic data samples di in d into the spectral components mi in m.
In general, problems like Equation (1) can be solved by a least squares approach, such as
m=(LHL)xe2x88x921LHd.xe2x80x83xe2x80x83(2) 
However, the matrix LHL may not be invertible for the Radon transform matrix L. Thus, Equation (1) is typically solved as a diagonally stabilized least squares problem given by
m=(LHL+kI)xe2x88x921LHd.xe2x80x83xe2x80x83(3) 
Here k is a pre-whitening factor or diagonal stabilization constant and I is the Nxc3x97N identity matrix. The diagonal stabilization matrix kI is added to the matrix LHL to avoid numerical instabilities in the matrix inversion. The superscript H indicates the Hermitian transpose of a matrix. The matrix (LHL+kI) to be inverted has Toeplitz structure and thus the inversion can be done using Levinson recursion.
Once the data has been transformed into the Radon domain as the model solution m, the primaries can be muted from m. Then, the multiples remaining in m can be transformed back into the spatial domain by the inverse transform
dm=Lmm,xe2x80x83xe2x80x83(4) 
where dm is the data vector containing substantially just multiples and mm is the model containing substantially just multiples. Now the multiples in multiple data vector dm can be subtracted from the data vector d to give substantially just the desired primaries. The matrix to be inverted, (LHL+kI), has Toeplitz structure and can be inverted using the Levinson recursion scheme, since the stabilization matrix kI is constant.
Hampson, D., 1986, xe2x80x9cInverse Velocity Stacking for Multiple Eliminationxe2x80x9d, J. Can. Soc. of Expl. Geophys., 22(1): 44-55, describes a parabolic Radon transform. It is used to attenuate multiples in a computer-efficient way. This is the method described above. However, it only works well with large move outs. These large move outs are, for example, twenty or more milliseconds at the farthest usable offset.
Sacchi, M. D. and Ulrych, T. J., 1995, xe2x80x9cHigh-Resolution Velocity Gathers and Offset Space Reconstructionxe2x80x9d, Geophysics, 60, 4, 1169-1177, describe a high-resolution Radon transform for attenuating multiples. The high-resolution transform constrains the Radon spectra to be sparse, using a sparseness prior for the curvature direction. This could also be called a re-weighted iterative approach. However, the high-resolution transform can show extra artifacts, which may compromise amplitude preservation. In addition, amplitudes may be affected if the induced sparseness is too strong.
Hunt, L., Cary, P., and Upham, W., 1996, xe2x80x9cThe Impact of an Improved Radon Transform on Multiple Attenuationxe2x80x9d, 66th Ann. Internat. Mtg.: Soc. of Expl. Geophys., 1535-1538, shows that a high-resolution parabolic Radon transform based on Sacchi and Ulrych""s 1995 publication improves the separation of primaries from multiples that have a small move-out difference. However, the techniques still suffers from the same deficits as Sacchi and Ulrych does.
Zwartjes, P. and Duijndam, A., 2000, xe2x80x9cOptimizing Reconstruction for Sparse Spatial Sampling, 70th Ann. Internat. Mtg.: Soc. of Expl. Geophys., 2162-2165, show that high-resolution Fourier regularization shows a better reconstruction in large gaps than standard least squares regularization.
In a high-resolution Radon transform, the diagonal matrix kI in the diagonally stabilized least squares approach given by Equation (3) is replaced by a Nxc3x97N diagonal matrix Q. The matrix Q is a stabilization matrix. Thus, Equation (3) is replaced by
m=(LHL+Q)xe2x88x921LHd.xe2x80x83xe2x80x83(5) 
The matrix Q in a high-resolution Radon transform is more than just a stabilization matrix. The matrix Q is now also used to constrain the solution m to the dominant curvatures q in the Radon spectral domain. Thus, the matrix Q will now be referred to more generally as a constraint matrix, rather than just a stabilization matrix.
Crucially, the constraint matrix Q is dependent upon the solution for the Radon domain m, and, conversely, that solution m is dependent upon Q. Thus, Equation (5) is a nonlinear problem and hence difficult, that is, expensive, to solve. Making the Radon transform approach to attenuating multiples more useful requires finding more computationally efficient methods for calculating or estimating the constraint matrix Q.
Sacchi and Porsani, 1999, xe2x80x9cFast high-resolution parabolic Radon transformxe2x80x9d, 69th Ann. Internat. Mtg.: Soc. of Expl. Geophys., 1477-1480, propose to solve this problem, such as given in Equation (10), as an iterative linear problem. The constraint matrix Q is iteratively derived using the solution for m from the previous step. FIG. 1 shows a flowchart illustrating the processing steps of the method of Sacchi and Porsani (1999) for processing seismic data.
At step 101, a seismic data set d is selected for seismic data processing. The seismic data set d is typically a gather of recorded seismic traces.
At step 102, the seismic data set d selected in step 101 is transformed from the (x, t) domain to the (x, xcfx89) domain. This transformation generates a transformed data set dF. The transformation defines a set of frequencies xcfx89 in dF. The transformation is typically done by a Fourier transform, preferably by a fast Fourier transform.
At step 103, a set of Radon spectral components are selected for a Radon transform. The spectral components are a set of curvatures q.
At step 104, a frequency co is selected from the set of frequencies obtained in the transformed data set dF in step 102. The frequencies would typically be selected in sequential order, starting at the lowest frequency.
At step 105, a matrix L representing a Radon transform into the set of curvatures q selected in step 103 is constructed. The matrix L satisfies Equation (1) for a Radon model solution m as defined in that discussion. Thus, the solution m is defined as a function of the curvatures q at each frequency xcfx89. At each frequency xcfx89, the solution m will be solved as an iterative series of linear problems giving mn, for n=0, 1, . . . , for however many iterations are required to arrive at a satisfactory solution.
At step 106, an initial Radon model solution mo is calculated for the frequency a) by the following straightforward approach,
m0=LHd.xe2x80x83xe2x80x83(6) 
At step 107, a constraint matrix Qi-1 for the Radon transform matrix L is derived from the Radon model solution mi-1 calculated at the previous iteration in either step 106 or 108. The matrix Qi-1 is defined as is the matrix Q in Equation (5), above.
At step 108, another iterative Radon model solution mi is calculated at the frequency xcfx89, using the value for the constraint matrix Qi-1 derived in step 107, by the following version of Equation (5),
mi=(LHL+Qixe2x88x921)xe2x88x921LHd.xe2x80x83xe2x80x83(7) 
At step 109, it is determined if another iteration is required for another iterative solution mi to obtain a satisfactory solution for m. Any suitable criteria may be used to determine if the current value of mi is satisfactory. If the answer is no, then the process returns to step 106 to begin another iteration. If the answer is yes, then the process continues to step 110.
At step 110, it is determined if all the frequencies co desired have been selected. If the answer is no, then the process returns to step 104 to select another frequency. If the answer is yes, then the process continues to step 111.
At step 111, the data is transformed from the (q, xcfx89) domain to the (q, tau) domain. Here tau is the zero offset intercept time. The transformation is preferably done by an inverse fast Fourier transform.
The method of Sacchi and Porsani""s 1999 publication uses an iterative linear scheme at each frequency to derive the constraint matrix Q from the solution of the model m from the previous iteration, but at the same frequency.
The international application, No. PCT/IB00/00659, published as International Publication Number WO 00/67045, xe2x80x9cA High resolution Radon Transform Seismic Traces Processingxe2x80x9d, by Compagnie Generale de Geophysique, proposes a non-iterative approach. They propose to derive the constraint matrix Q from the solution m at the previous frequency. They describe a constrained high-resolution parabolic Radon transform method. The method involves performing a constrained high-resolution Radon decomposition at various frequencies in seismic traces, wherein the Radon decomposition at a given frequency is constrained as a function of the Radon decomposition at at least a lower frequency. In particular, it describes successively performing a Radon decomposition at various sparse frequencies, wherein the Radon decomposition at a given frequency is constrained as a function of the Radon decomposition at the previous frequency. This has the computational advantage of avoiding the re-weighted iterative approach of Sacchi and Ulrych (1995) or Sacchi and Porsani (1999). However, the constraint matrix Q used in this Radon transform is still frequency dependent. FIG. 2 shows a flowchart illustrating the processing steps of a method of the Compagnie Generale de Geophysique patent application publication for processing seismic data.
At step 201, a seismic data set d is selected for seismic data processing. The seismic data set d is typically a gather of recorded seismic traces.
At step 202, the seismic data set d selected in step 201 is transformed from the (x, t) domain to the (x, xcfx89) domain. This transformation defines a set of frequencies xcfx89. The transformation is preferably done by a fast Fourier transform.
At step 203, a set of Radon spectral components are selected for a Radon transform. The spectral components are a set of curvatures qi for i=xe2x88x921, 0, 1, . . . , N.
At step 204, a frequency Co is selected from the frequencies obtained in the seismic data set d from its transformation in step 102. The frequencies would typically be selected in sequential order, starting at the lowest frequency.
At step 205, a matrix L representing a Radon transform into the set of curvatures qi selected in step 203 is constructed. The matrix L satisfies Equation (1) for a Radon model solution m as defined in that discussion. Thus, the solution m is defined as a function of the curvatures qi for each frequency xcfx89i. At each frequency xcfx89i, the solution m(xcfx89i) will be solved as a linear problem.
At step 206, an initial Radon model solution m(xcfx89xe2x88x921) should be calculated for the frequency xcfx89xe2x88x921. However, the Compagnie Generale de Geophysique patent application publication does not mention how they would do this. Therefore, it can only be assumed that an initial solution m(xcfx89xe2x88x921) is somehow obtained. One simple, straightforward approach would be the following equation,
m(xcfx89xe2x88x921)=LHd(xcfx890).xe2x80x83xe2x80x83(8) 
At step 207, a frequency xcfx89i is selected from the set of frequencies selected in step 204. The frequencies are chosen in sequential order, starting at the lowest.
At step 208, a constraint matrix Qixe2x88x921 for the Radon transform matrix R is derived from the Radon model solution m(xcfx891xe2x88x921) calculated at the previous frequency xcfx89ixe2x88x921 in either step 206 or 209. The matrix Qixe2x88x921 is defined as is the matrix Q in Equation (5), above.
At step 209, a Radon model solution m(xcfx89) is calculated at the frequency xcfx89i, using the value for the diagonal constraint matrix Qixe2x88x921 derived at the previous frequency xcfx89ixe2x88x921 in step 208, by the following version of Equation (10),
m(xcfx89i)=(LHL+Q(xcfx89ixe2x88x921))xe2x88x921LHd.xe2x80x83xe2x80x83(9) 
At step 210, it is determined if all the frequencies have been selected. If the answer is no, then the process returns to step 204 to select the next frequency. If the answer is yes, then the process continues to step 211.
At step 211, the data is transformed from the (q, c) domain to the (q, tau) domain. Here tau is the zero offset intercept time. The transformation is preferably done by an inverse fast Fourier transform.
The method of the Compagnie Generale de Geophysique, described in FIG. 2, is non-iterative and thus more computationally efficient than the method of Sacchi and Porsani (1999), described in FIG. 1. In addition, by using a lower frequency as a constraint or stabilization on the next frequency, this method is better at handling aliased data. However, the matrix Q is frequency dependent. The approximation for Q is calculated for each frequency in the gather.
Hargreaves, N. and Cooper, N., 2001, xe2x80x9cHigh-Resolution Radon Demultiplexe2x80x9d, 71st Ann. Internat. Mtg.: Soc. of Expl. Geophys., 1325-1328, show that the high-resolution transform can cause smearing of energy in the transform domain in the direction in which no sparseness is imposed, i.e. the time direction. This is seen in their FIG. 2. It was also apparent in FIG. 4a of Sacchi and Ulrych""s 1995 publication, discussed above. This smearing raises the question of amplitude preservation. In addition, amplitude preservation may not be optimal if the induced sparseness would be too strong. For example, if a perfectly parabolic event with AVO variation is focused onto a single curvature parameter, then the AVO variation is lost.
Therefore, a need exists in the field of seismic data processing for a high-resolution Radon transform method for attenuating multiples that is potentially more computationally efficient. Thus, a need exists for a method to calculate a constraint matrix that is both non-iterative and frequency independent. Additionally, the method should preserve amplitudes, yet avoid smearing energy in the time direction.
The invention is a method for processing a seismic data set. A seismic data set is processed by selecting a set of gathers in the seismic data set. A frequency independent constraint matrix is calculated in a first gather. A model solution in at least one second gather is calculated by using the constraint matrix in a high-resolution Radon transform.
In an alternative embodiment, the primaries are attenuated in the model solution, generating a multiples model solution. The multiples model solution is inverse transformed from the Radon domain back to the spatial domain, generating a multiples data set. The multiples data set is subtracted from the seismic data set, generating a substantially multiple-free data set of primaries.